Quenches in a quasi-disordered integrable lattice system: 
Dynamics and statistical description of observables after relaxation 
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We study the dynamics and the resulting state after relaxation in a quasi-disordered integrable lattice system 
after a sudden quench. Specifically, we consider hard-core bosons in an isolated one-dimensional geometry 
in the presence of a quasi-periodic potential whose strength is abruptly changed to take the system out of 
equilibrium. In the delocalized regime, we find that the relaxation dynamics of one-body observables, such as 
the density, the momentum distribution function, and the occupation of the natural orbitals, follow, to a good 
approximation, power laws. In that regime, we also show that the observables after relaxation can be described 
by the generalized Gibbs ensemble, while such a description fails for the momentum distribution and the natural 
orbital occupations in the presence of localization. At the critical point, the relaxation dynamics is found to be 
slower than in the delocalized phase. 
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I. INTRODUCTION 



The nonequilibrium dynamics of isolated integrable quan- 
tum systems is constrained by a large number of conserved 
quantities, which generally preclude relaxation to thermal 
equilibrium lfll-l23ll. This may affect current experiments that 
are realized in one-dimensional (ID) and quasi- ID geome- 
tries close to integrable points [12411 and future technological 
devices. As such, this phenomenon cannot be considered as 
purely academic anymore. Advances in controlling and ma- 
nipulating highly isolated quantum gases in low dimensions, 
and at very low temperatures, have made it possible to study 
in great detail the relaxation dynamics following an abrupt 
change of some of the system's parameters ||25[ l26ll . so that 
questions related to the lack of thermalization can now also be 
addressed experimentally. For example, in Ref. ll25ll . it was 
experimentally shown that the relaxation dynamics of one- 
dimensional atomic Bose gases does not necessarily lead to 
a thermal momentum distribution of the atoms. 



Soon after the experimental finding in Ref. 025II . it was 
shown in Ref. JH] that expectation values of few-body ob- 
servables in isolated integrable systems after relaxation can be 
predicted by generalized Gibbs ensembles (GGEs). GGEs are 
constructed by maximizing the entropy i27[|28ll . while satisfy- 
ing constraints imposed by the constants of motion that make 
the system integrable. Interestingly, the mechanism that leads 
to thermalization in non-integrable systems, namely, eigen- 
state thermalization HH-IUlL can be generalized to the inte- 
grable case in the sense that most eigenstates that are close 
not only in energy but also in their distribution of conserved 
quantities share the same expectation values of few-body ob- 
servables JH. This allows one to understand why the GGE 
description works. Its validity after relaxation has been tested 
in many different integrable quantum models |[ll-l23ll. and 
has been argued to be adequate for predicting prethermalized 
expectation values of observables IB3K35I1 in non-integrable 



quantum systems close to an integrable point 13( 

In relation to current ultracold-gas experiments (and to low- 
dimensional mesoscopic devices), one question that needs to 
be addressed is the fate of the GGE description when trans- 
lational invariance is absent in the system. Numerical calcu- 
lations for hard-core bosons in a box JU-H]] and in the pres- 
ence of a harmonic confining potential (relevant to optical 
lattice setups) |0,[l7t], have shown that the GGE indeed de- 
scribes observables after relaxation. However, a recent study 
of quenches in the quantum Ising chain has put forward the 
notion that "as soon as the translational invariance is broken, 
the GGE fails to apply" ll37ll . This was supported by calcula- 
tions of equal-time correlations after a quench in the presence 
of disorder. Since the general statement made in Ref. ll37ll 
is in contradiction with previous results lfll-l3l especially 
with those in the presence of a confining potential |Hl[l7l], nere 
we reconsider the question of whether the GGE description is 
valid in the absence of translational invariance. 

One important difference between the systems studied in 
Ref. ll37ll and those studied in Refs. OJ-Hl] is the inclusion of 
disorder in the former. Even in the presence of interactions, 
disorder can lead to localization Il38l - l4lll . and, in noninte- 
grable systems, localization can lead to lack of thermalization 
after relaxation following a quantum quench PH l42ll . The lat- 
ter can be understood to follow from the failure of eigenstate 
thermalization in the localized regime filTl . It is then natu- 
ral to expect that, in integrable systems, localization, and not 
necessarily the breaking of translational symmetry, may lead 
to a failure of the GGE description. This would follow from a 
failure of the generalized eigenstate thermalization . 

In order to separate the effects of breaking translational 
symmetry and localization in an integrable system, we study 
hard-core bosons in an incommensurate superlattice. This 
model exhibits a transition between an extended and a lo- 
calized phase at a finite strength of the superlattice potential 
ll43ll . and is to be contrasted with the case of uniform random 
disorder where localization occurs for any nonzero disorder 
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strength rf24ll . We show that in the extended phase, the GGE 
provides a correct description of one-body observables after 
relaxation, despite the lack of translational invariance. On the 
other hand, in the localized phase, the GGE is found to fail. At 
the critical point, a slower relaxation dynamics is seen to pre- 
clude the observation of stationary values of the observables 
for the largest system sizes. However, as long as the station- 
ary value is reached, the GGE provides a good description of 
observables after relaxation at the critical point. 

The exposition is organized as follows. In the next sec- 
tion (Sec. HD, we introduce the model and observables to be 
studied in the remainder of the paper. We also briefly discuss 
the computational approach utilized in our study, as well as 
the ensembles that are used to compare with the results after 
relaxation. In Sec. Hill we study the relaxation dynamics fol- 
lowing a sudden quench in the different regimes of the model. 
Section|lV]is devoted to a comparison of observables after re- 
laxation with the predictions of statistical ensembles, as well 
as a finite-size-scaling analysis that allows us to gain insight 
into the behavior in the thermodynamic limit. We also make 
contact with the results in Ref. 13711 by studying the behavior 
of off-diagonal one-particle correlations. The conclusions are 
presented in Sec. [V] 

II. MODEL, OBSERVABLES, AND ENSEMBLES 

Our study is performed within the Aubry-Andre model [R3I1 
for hard-core bosons in a one-dimensional lattice with open 
boundary conditions. The Hamiltonian reads 

L-l L 

# = - f E (£py+i+H.c.)+A£cos(27tC7j + (p)^, (1) 

where the operator Pj (bj) creates (annihilates) a hard-core 
boson at site j, and hj = b^bj is the on-site occupation number 
operator, bj and bj obey the usual bosonic commutation rela- 
tions, i.e., j] = 8ij, but satisfy a constraint bj = b' 2 = 0, 
which forbids multiple occupancy of the lattice sites. The hop- 
ping parameter is denoted by t (we set t = 1, h = 1 throughout 
this work), L is the number of sites, and we consider only sys- 
tems in which the number of particles (AO is N = L/2 (half 
filling). By selecting a to be an irrational number, we gener- 
ate a quasiperiodic potential whose strength is controlled by 
the parameter A. In our study, we choose a to be the inverse 
golden ratio, a = (V5 — l)/2, a choice motivated by the fact 
that the golden mean is considered to be the most irrational 
number 14411 . (p allows the phase of the potential to be shifted, 
and will be used later to average over different realizations in 
our finite systems. For most of our work, we set <p = 0. 

Despite the quadratic form of Eq. ([T}, it cannot be di- 
rectly diagonalized because of the on-site constraints forbid- 
ding multiple occupancy of the lattice sites. This can, how- 
ever, be circumvented by mapping the ID hard-core boson 
Hamiltonian onto a spin- 1/2 chain via the Holstein-Primakoff 
transformation |45|], and then map ping the spin- 1/2 chain onto 
noninteracting spinless fermions Il46ll via the Jordan- Wigner 



transformation firTl . The resulting Hamiltonian maintains the 
form in Eq. (HJ but with the hard-core operators replaced by 
fermionic ones. It then follows that the spectrum, as well as 
thermodynamic and density-related properties, are the same 
for hard-core bosons and non-interacting spinless fermions. 

The Aubry-Andre model ll43ll is known to undergo a local- 
ization transition at a critical X c = 2t. For A < A c , all single- 
particle states are extended, i.e., Bloch-like, states. Above 
the critical point, single-particle states are exponentially lo- 
calized with localization length % = ln(A) -1 114311 - Because of 
the mapping above, the same holds true for hard-core bosons. 
This implies that the ground state of the latter undergoes a 
superfluid-insulating transition as X c = 2f is crossed. In the 
localized phase, the ground state is a Bose glass lf24ll . 

In connection to optical lattice experiments, such as the 
ones carried out in Refs. ll25[|26ll . we are interested in study- 
ing two different one-body observables: the on-site density 
iij = (77,), and the momentum distribution function m^. is 
the diagonal part of the Fourier transform of the one-particle 
density matrix p !; - = (bjb •}, 

ij=l 

Additional information on the coherence properties of the sys- 
tem can be gained through the study of the natural orbitals 
0" and their occupations 7] a , defined through the eigenvalue 
equation 

tpy07 = %0f. (3) 

7=1 

In homogeneous periodic systems, the natural orbitals are 
plane waves and their occupations coincide with the momen- 
tum distribution function, so and rj a give the same infor- 
mation about the system. However, once translational invari- 
ance is broken these two quantities become different. Out 
of equilibrium, they can even give apparently inconsistent re- 
sults. For example, during the expansion of a hard-core boson 
gas its momentum distribution function becomes identical to 
that of noninteracting fermions, which may be taken as an in- 
dication that the system lacks coherence 14811 . However, the 
occupation of the natural orbitals is very different from the one 
of fermions; many orbitals remain highly populated, which re- 
veals the bosonic character of the out-of-equilibrium gas ll48ll . 
In addition, in higher-dimensional interacting systems, if the 
occupation of the highest occupied natural orbital scales with 
the total number of particles, then one can say that the system 
exhibits Bose-Einstein condensation [49, 50J. 

In equilibrium, the properties of hard-core bosons, modeled 
by Eq. ([T]), have been studied in detail in the ground state f5ll 
l52ll and at finite temperature ll53ll . Here, our goal is to examine 
the dynamics after the system is taken out of equilibrium by a 
sudden change of A (A/ — > Xf). The initial state ^i) is taken 
to be the ground state of Hj [Eq. ([TJ with A = A/] and the 
evolution is studied under Hp [Eq. ([Tji with A = Xf\: 

|V(T))=e- i * t |Y / ). (4) 
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To study the time evolution of the observables introduced 
above, we follow a computational method based on the Bose- 
Fermi mapping and the use of properties of Slater determi- 
nants. This method has been explained in detail in Ref. fsill . 
so we do not reproduce it here. It allows one to calculate 
each matrix element p,y (at any given time x) in terms of 
the determinant of an (N + 1 ) x (N + 1 ) matrix, which results 
from the product of two matrices with sizes (N +1) x L and 
L x (N + 1). The computation time of the entire one -particle 
density matrix essentially scales as L 2 (N + l) 3 (the matrix 
multiplication need not be done for every entry), which al- 
lows us to efficiently study the dynamics of systems of up to 
1000 lattice sites. 

We then contrast the time-averaged results for the observ- 
ables after relaxation with the predictions of statistical me- 
chanics. While the most relevant traditional ensemble to 
compare with would be the microcanonical one (because the 
time-evolving system is isolated), we instead use the grand- 
canonical ensemble (GE). This is because calculations in the 
former scale exponentially with system size, while, in the lat- 
ter, they scale as power laws. Within the GE, we can study 
very large lattices, in which we expect a good agreement 
between the predictions from different statistical ensembles 
ll55ll . The density matrix in the GE reads 



1 

Pge = -—exp 



k B T 



(5) 



where kg is the Boltzmann constant, N is the total number 
operator, and Zqe is the partition function, 



Z G e = Tr 



exp 



H-fiN 



(6) 



In order to compare the grand-canonical predictions for the 
observables to those obtained following the quantum dynam- 
ics, T and jU need to be chosen so that Tr[pGE#f] = E and 
Tr[p GE A>] = N, where E = ( l F/|H F | l F/) is the energy of the 
time-evolving system after the quench, which is conserved. 

In integrable hard-core-boson systems, in the absence of 
disorder or quasi-disorder, the grand-canonical (HJ-IH and mi- 
crocanonical Ql descriptions have been shown to fail to pre- 
dict the outcome of the relaxation dynamics for few-body ob- 
servables. Instead, the GGE has been proposed to be the ade- 
quate ensemble to deal with this problem JU]. The GGE den- 
sity matrix can be written as 



Pgge 



1 



Zgge 



exp 



" X m I f) 



(7) 



where I m are the conserved quantities, X m are their corre- 
sponding Lagrange multipliers, and Zqge is the partition func- 
tion, 



Zgge — Tr 



exp 



(8) 



The Lagrange multipliers need to be selected so that the ex- 
pectation values of the conserved quantities in the GGE are the 



same as in the initial state, i.e., Tr[pGGE^n] = (^/l^ml 1 ?/)- For 
hard-core bosons, where the conserved quantities are taken 
to be the projection operators to the single-particle eigenstates 
of the fermionic Hamiltonian to which Eq. ([T} can be mapped, 
the Lagrange multipliers can be written as yj 



Xm = In 



I-OP/I/mIT/) 

OF/l/ml*/} 



(9) 



In order to calculate the expectation value of the one- 
particle density matrix in the grand-canonical ensemble, 



,GE 



Tr 



b;bjP GE 



, and in the GGE, p 



GGE 



Tr 



bj bjPGGE 



(note that the GGE is also grand-canonical), we use the ap- 
proach introduced in Ref. Il55ll . The grand-canonical calcu- 
lations, similarly to the ones carried out for studying the dy- 
namics, use the Bose-Fermi mapping and properties of Slater 
determinants. The computation time of the entire one-particle 
density matrix in this case scales as L 5 f55ll . 



III. TIME EVOLUTION 

To probe the relaxation dynamics after the quench, we cal- 
culate the normalized difference 80(z) (where O stands for n, 
m, rj) between the expectation value of observables at differ- 
ent times and their long-time average Oy. 80(x) is defined 
as 



80{x) 



E; 0j(x)-O? 



(10) 



[Note that j is a dummy variable that stands for i (in n,), k 
(in rrik), and a (in rj a )]- If observables relax to stationary 
values, 80(x) will fluctuate about a time-independent value. 
This value, as well as the amplitude of the time fluctuations 
about it, is expected to be finite for finite systems but should 
vanish in the thermodynamic limit. We note that OJ is taken 
to be an average over a variable size time interval that contains 
the longest times that we have simulated. In the event that the 
observable has not relaxed by then, 80(x) will make it evident 
as it will not become stationary. 

In Fig.Q] we show results for 80(t) in a set of quenches in 
which the initial state is the ground state of Eq. dTJ with A/ = 
(i.e., a superfluid state) and Xf is below (Xf = 1), at (Xf = 2), 
and above (Xf = 3,4) the localization transition. Results are 
presented for three different system sizes (L = 10, 100, and 
1000, from top to bottom in each panel). In Figs. 02 a )-02 c X 
one can see that all three observables in the quench terminat- 
ing in the extended phase exhibit a clear relaxation dynamics 
in which 80(x) decreases as time passes, and then fluctuates 
about a finite time-independent value. Both the finite time- 
independent value and the amplitude of the fluctuations are 
seen to decrease with increasing system size. 

The quench towards the critical point (Xf = 2) [Figs.QJd)- 
|TJf)] exhibits a different dynamics. As the system size in- 
creases beyond 100 sites, the three observables considered 
here do not reach a clear stationary value during the times 
studied (up to x = 10 6 for «, and x = 5.37 x 10 4 for and 
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8n 5m 8r\ 




FIG. 1: (Color online) Relaxation dynamics of n,-, m^, and r\ a as they approach the long-time average in a quench A/ = — s> Xp = 1 (a)-(c), 
2/ = -> A F = 2 (d)-(f). A/ = — > = 3 (g)-(i), and A/ = A F = 4 (j)-(l), for systems with 10, 100, and 1000 lattice sites (from top 
to bottom in each panel). The time averages are computed as follows: Since »; is computationally less expensive than and rj a , for that 
observable we simulated longer times and averaged over 9000 steps with X <= [10 5 , 10 6 ] for all lattice sizes. For and rj a , we averaged over 
900 steps with z £ [10 4 , 10 5 ] for L = 10 and L = 100, and over 437 steps for L = 1000 with T e [10 4 , 5.37 x 10 4 ] (in the plots, T max = 5.37 x 10 4 ). 



r\ a ). This can be understood as the critical point is known to 
be very special. The single-particle spectrum becomes a Can- 
tor set (the bands acquire zero measure), and the gaps form 
a devil's staircase IprJ . Such a peculiar spectrum seems to 
render dephasing ineffective in these systems. Our finding 
implies that, at the critical point, stationary values of the ob- 
servables may be more difficult to observe experimentally. 

Finally, the quench towards the localized regime 
[Figs. mg)-Ei) for X F = 3 and Figs. [Uj)4Hl) for X F = 4] 
does lead to stationary values for and rj a . Note that 



and rja exhibit dynamics that are qualitatively similar to the 
one observed in the quench A/ = — > Aj? = 1, namely, the 
stationary values of 8m and 8rj (and the fluctuations about 
them) decrease with increasing system size. «,, on the other 
hand, exhibits a different behavior. Because of localization 
in real space, 8n becomes lattice-size independent, i.e., it 
remains finite in the thermodynamic limit. In that case, the 
only effect that increasing L has is to reduce the amplitude of 
the time fluctuations of 8n about the stationary value. 

We have also studied quenches starting from different ini- 
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FIG. 2: (Color online) As Fig.Q]but for quenches from Xj = 8, i.e., from deep inside the Bose-glass phase. 



tial states that are eigenstates of Eq. ((TJ, and even from 
the ground state of commensurate superlattices such as the 
ones studied in Refs. HKjjUlj]], finding a qualitatively sim- 
ilar dynamics to the one depicted in Fig. Q] As an exam- 
ple of a different initial state, in Fig. [2] we report results in 
which the quenches start from the ground state of Hamilto- 
nian (Q~|i deep inside the Bose-glass phase (A/ = 8). Figure 
[2] shows that the dynamics is indeed very similar to that re- 
ported in Fig. Q] The only apparent difference is that for 
quenches within the Bose-glass phase (A/ = 8 — > Xp = 3 and 
A/ = 8 — > Xp = 4), the stationary value of 8n is smaller than 
in the quenches from the superfiuid phase to the Bose-glass 
phase (A/ = —¥ Xp = 3 and A/ = — > Xp = 4). For the for- 
mer, we find <5n 8 ^ 3 H w 0.06 and <5n 8 ^ 4 H « 0.04 while 



for the latter <5«°^ 3 (°°) 5n°^ 4 (oc) « 0.15. This is under- 
standable as 8n(0) is already smaller in quenches starting in 
the Bose-glass phase than in those starting in the superfiuid 
phase. 



Approach to the stationary values 

In a recent numerical study of the relaxation dynamics of a 
disordered nonintegrable fermionic system with short-range 
interactions and random long-range hopping, it was found 
that, in the extended phase, observables exhibit a power-law 
approach to their thermal expectation values ll4Hl , Power-law- 
like relaxation dynamics was also seen in recent optical lattice 
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FIG. 3: (Color online) (a) 5m* vs X for <p = 0, as well as after averag- 
ing over 1000 random values of <p (uniformly distributed in [0,27f]), 
in systems with 100 lattice sites. The fits to power-law and expo- 
nential behavior were done over the interval T 6 [1,40] (a vertical 
line marks T = 40), which contains 1200 data points, (b) Sm^ vs T 
for <p = in a system with 1000 lattice sites. The fits to power-law 
and exponential behavior were done over the interval T € [1,600] (a 
vertical line marks T = 600), which contains 230 data points. 



experiments with a clean system in a one-dimensional geom- 
etry l26ll . These results are to be contrasted with the expo- 
nential approach expected in generic nonintegrable systems. 
Since both studies I26ll4lll were limited to small lattice sizes, 
and no extensive scaling analysis could be performed, it is not 
clear how these findings are affected by finite-size effects. 

The dynamics depicted in Figs. Q] and [2] for three system 
sizes, which are a decade away from each other, provides a 
clearer picture of the role of finite-size effects. We indeed find 
indications of power-law relaxation, as it is apparent in the 
plots that the time interval over which a power-law-like behav- 
ior is seen increases with system size. We explicitly show this 
in Fig. [3] where we compare the relaxation process for sys- 
tems with 100 and 1000 lattice sites. In the former [Fig. [2 a)], 
both power-law and exponential decay provide a reasonable fit 
to the data. In the latter [Fig.|3jb)], where power-law behav- 
ior is apparent for about three decades, a fit to an exponential 
decay is clearly inconsistent with the data. Hence, our results 
provide another example of a system in which, whenever re- 
laxation takes place, the relaxation dynamics follows a power 
law. To what extent power-law-like relaxation is generic to 
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FIG. 4: (Color online) Finite-size scaling of 5«(°°), 8m(<x>), and 
<5?7(°°) for the quenches studied in Figs. [7] and [2] The dashed 
lines are power-law fits leading to 8n(°°) oc L~ 0,49 , 8m(°°) oc L -052 , 
and 5rj(°°) oc L~ a51 in (a), gn(°°) oc ZT - 50 , 5m(°°) °c L~ 052 , and 
577 H °= ^~ a51 in (b), 8n{°o) oc L~ - 25 in (c), Sn(°°) « L~ 0,26 
in (d), 5»(°°) oc IT - 01 , Sm(°°) oc L~ 0,43 , and 577(00) oc L~ 0,49 in 
(e), 8n(°o) oc L" 01 , SmM oc L -0 ' 51 , and 577(00) oc L~ 0,48 in (f), 
5n(oo) oc L°, 5m(oo) oc L M , and 577(00) oc L~ - 49 in (g), and 
5;i(oo) oc L°, 5m(oo) oc L-° 50 , and 577(00) oc L~ 0,48 in (h). The 
power-law fits were done using the data for systems with between 
100 and 1000 lattice sites (11 data points). 



the dynamics of isolated quantum systems, especially nonin- 
tegrable ones, is a topic that deserves further attention. 

Since we are dealing with finite lattice sizes with open 
boundary conditions, we have also studied the effect that aver- 
aging over different phases <p [see Eq. (Q]i] has on our results. 
A typical outcome of such an average is depicted in Fig. |3a), 
for 1000 different values of (p distributed uniformly in [0,2n]. 
The average over different phases can be seen to reduce time 
fluctuations after relaxation, but leaves the results for the ap- 
proach to the stationary value almost unaffected. 

Another important question to be answered, which is of 
special interest to current experiments with ultracold gases, 
is how long it takes for observables to reach the stationary 
values. Given the strong indications found above that the re- 
laxation dynamics follows a power law, the times at which sta- 
tionary values are attained will be determined by how 5n(°°) ; 
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8m(°°), and 8r) (oo) (here, should be understood as a long 
time after relaxation) scale with system size. In Fig. @] we 
show the scaling of those quantities in the quenches analyzed 
in Figs. [T]and E] Figures Q3a), G3b), andHte)-QIf) show that, 
away from the critical point, the scaling of 8m(°°) and 8r](°o) 
is close to 1/VX, and a similar scaling is seen for 8n(°°) in 
quenches to the extended phase [Figs. [TJa) and [Tib)]. Such 
a scaling has been proven to provide a bound for the normal- 
ized time variance of observables that are quadratic in Fermi 
operators in noninteracting-fermion models |[57ll . but we find 
it to be also applicable to more general observables in inte- 
grable systems. As discussed before, in quenches to the lo- 
calized regime, 8n(°°) becomes independent of system size. 
Also, the slow relaxation dynamics of m and rj at the critical 
point precludes the observation of a clear scaling for 8m(°°) 
and 8r](°°) [Figs. [He) and[TJd)], while the scaling of is 
close to 1 /L 1 / 4 . The scalings of 5n(°°) at the critical point and 
in the localized regime violate the bound proven in Ref. ll57ll . 

A power-law approach of 5«(t), <5to(t), and 5tj(t) to the 
stationary values, together with a power-law scaling of 8n(°°), 
8m(°°), and 5rj(°°) with system size, implies that the time at 
which stationary values are attained increases as a power law 



with system size. This means that measuring densities and 
momentum distribution functions in experiments is advanta- 
geous with respect to directly measuring two-point correlation 
functions. After relaxation, the values of the latter have been 
shown to be exponentially small compared with the distance 
between the points Q I22D and, as such, the time it takes for 
those correlations to relax to the stationary values increases 
exponentially with the distance between the points 12311 . 



IV. DESCRIPTION AFTER RELAXATION 

After discussing the relaxation dynamics, we focus on the 
description of the observables after relaxation. In generic 
(non-integrable) quantum systems, one expects the dynam- 
ics to lead to thermalization, namely, to expectation values 
of observables that are equal to those of a system in thermal 
equilibrium. Because of thermodynamic universality, this is 
expected to be true whenever the isolated system and its ther- 
mal equilibrium counterpart share the same mean energy and 
number of particles 029143211 . independently of the initial state 
in the former. 
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FIG. 5: (Color online) Density in the first ten sites (a),(d), momentum distribution function (b),(e), and natural orbital occupations (c),(f) for 
quenches in which the initial state is the superfluid ground state of a system with Xj = while Xf = 1 (a)-(c), Xf = 4 (d)-(f), and for L = 1000. 
We present results for the observables in the initial state, the long-time average [calculated between t = 10 5 and T = 10 6 for n, (9000 steps), 
and between T = 10 4 and T = 5.37 x 10 4 (437 steps) for mj and r/ a ; see the caption of Fig.fJJ, as well as within the GE and the GGE. Note 
that, except 8n(°°) for Xf = 4, 5n(<»), 5m(°°), and 5t](°°) are very small for L = 1000 (see Fig.0. In addition, we have checked that all time 
averages are well converged. 



In Fig. [5j we show results for nu m^, and r\ a for quenches 
from initial states with A/ = 0, and Xf = 1 and 4. For all 
quantities, we report their values in the initial state, the long- 
time averages, and within the GE and the GGE. The plots for 
the density in the initial state [Figs. [3 a) and [3d)] make ev- 
ident that, despite the presence of open boundary conditions, 



at T = the density is constant throughout the system. This is 
because of the particle-hole symmetry of the model. After the 
quench, this is not true anymore and the density becomes time 
dependent and inhomogeneous. The time-averaged result for 
the density after relaxation and the predictions of the GGE 
are indistinguishable from each other for Xf = 1 in Fig. [3 a) 
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FIG. 6: (Color online) As Fig.|5]but for quenches from Aj = 8, i.e., from deep inside the Bose-glass phase. 



and Xf = 4 in Fig. |5jd). The predictions of the GE are dif- 
ferent from the outcome of the relaxation dynamics in both 
Two other identifying properties of the initial state, which 
signal the existence of off-diagonal quasi-long range corre- 
lations, are the presence of a sharp peak in at k = 
[Figs.Eb) andHe)] and in r\ a at a = [Figs. Etc) andEtf)]. 
The quenches can be seen to lead to a dramatic decrease of 
the height of those peaks after relaxation, which is similar to 
the effect of finite temperature in equilibrium systems ll55ll . 
For nik and rj a , a stark contrast can be observed between the 
results obtained for the quench A/ = — >• Xf = 1 and those ob- 
tained for the quench A/ = — > Xf = 4. While, in the former, 
the time-averaged results and the GGE predictions are almost 
indistinguishable from each other, the same is not true for the 
latter. This suggests that the transition to localization plays an 
important role in the description after relaxation. In addition, 
the thermal values for both observables in the GE are clearly 
different from the results after relaxation. 

Qualitatively, we have obtained a very similar picture to the 
one gained through Fig. El for what happens after relaxation 
in the extended and localized regimes, for a wide range of dif- 
ferent initial states. Among those, we considered ground and 
excited states of hard-core-boson Hamiltonians in the form of 
Eq. ([TJ but with different local potentials, including period-2 
superlattices I0 [ill |2j]]. In Fig. El we show results for the 
case in which the initial state is the ground state of Eq. (Q3 
with A = 8. In contrast to the case with A/ = 0, for A/ = 8 
the initial state is deep inside the Bose-glass phase where the 
density is inhomogeneous [Figs. Eta) andEtd)] and the sys- 
tem lacks coherence. The latter is reflected by the almost flat 
initial momentum distribution [Figs.Etb) andEt e )L Localiza- 
tion in this regime is revealed by the natural orbital occupa- 
tions [Figs.Et c ) andEtf)]; which is nearly 1 for the first 500 



quenches. 

orbitals (there are 500 particles in the system), i.e., the bosons 
in this many-body system can be seen as single particles lo- 
calized within a few sites. This picture is confirmed by the 
form of the natural orbital wave functions (not shown). 

After the relaxation dynamics following the quenches A/ = 
8 — > Xp = 1 and A/ = 8 — > Xp = 4, one can infer from Fig. El 
[panels (b), (c), (e), and (f)] that one-particle correlations are 
enhanced from the ones in the initial state. This follows as 
the height of the zero-momentum occupations increases, the 
zero-momentum peaks become narrower, and the occupation 
of the lowest natural orbitals depart from 1 . This is very differ- 
ent from what happens in the quenches A/ = — > Xf 7^ de- 
picted in Fig. El where one-particle correlations are reduced. 
Despite this contrast, we find that the GGE results are al- 
most indistinguishable from the time-averaged ones for all ob- 
servables in quenches A/ = 8 — > Xp = 1, while for quenches 
A/ = 8 — > Xp = 4 only the density and are accurately de- 
scribed by the GGE. In the latter quench, the GGE fails to 
describe the natural orbital occupations, pointing once again 
towards the role of localization. 



Scaling with system size 

Even more important than the actual differences seen in 
Figs.E]andE]between the long-time averages and the predic- 
tions of statistical ensembles (GE and GGE) is how those dif- 
ferences scale with increasing system size (L = 1000 in those 
figures). One could imagine, for example, that while the dif- 
ferences between the time averages and the GE are large for 
finite systems they may disappear in the thermodynamic limit. 
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FIG. 7: (Color online) Finite-size scaling of A/?i GGE ' GE ) and 
Arj GGE(GE) for the quenches studied in Figs. Q] and [2] The dashed 



lines in (a)-(e) are power-law fits leading to Am' 



GGE 



At] 



in (a), and Am' 
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and Arj 
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-0.99 



and 
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in (b). Up to 100 sites, the time average was taken over 900 steps 
with T 6 [10 4 , 10 5 ]. For all other system sizes, the time average was 
taken over 437 steps with x e [10 4 ,5.37 x 10 4 ]. 



Another possibility is that the differences between the time 
averages and the GGE are small for the quenches and sys- 
tem sizes studied here but they may not vanish in the ther- 
modynamic limit, which would invalidate the GGE descrip- 
tion for thermodynamic systems. Cases in which integrable 
systems seemed to behave thermally, but failed to exhibit the 
required scaling with system size, were recently studied in 
Refs. CHH. 

In order to study the scaling of the discrepancies between 
the time averages and the statistical predictions, we compute 
the normalized differences AO between the long-time average 
of the observables Of and the ensemble predictions q gge< - ge ^ 



A0 GGE(GE) = 



Of - o 



GGE(GE) 

j 



(ID 



Note that O stands for n, m, and rj, and j is a dummy vari- 
able that stands for ;', k, and a, respectively. This quantity is 
defined in the same spirit as 80 in Eq. (fTOb . 

In Fig. El we show the scaling of Ara GGE(GE) and At7 gge(ge) 



for the quenches studied in Figs. Q] and [2] Apparent differ- 
ences can be seen between the scalings when Xf lies in the ex- 
tended, critical, and localized regimes. Different initial states, 
on the other hand, lead to qualitatively similar behavior of 
A(9 GGE ' GE \ i.e., Xf is the parameter that determines how the 
outcome of the relaxation dynamics compares to the predic- 
tions of statistical ensembles. 

In quenches terminating in the extended phase [Xf = 1, 
Figs. EI a) andElb)], one can see that Aw GGE and At] gge ex- 
hibit a power-law decrease with increasing system size. The 
small oscillations in Am GGE , seen in Fig. Elb) for the largest 
system sizes, are due to the small values of this quantity. They 
depend on the exact time intervals and number of time steps 
used in the time averages. Hence, such oscillations are an ar- 
tifact of our numerical calculations and are not expected to be 
present if one takes the infinite time averages used in previous 
works iHH]]], which are not available here. Am GE and Arj GE , 
on the other hand, exhibit a clear saturation to finite values 
with increasing system size. From these scalings, we con- 
clude that the GGE correctly describes and r\ a after relax- 
ation, despite the absence of translational invariance and the 
presence of disorder. On the contrary, the GE fails to describe 
those observables, which makes evident that these systems do 
not thermalize in the traditional sense. 

Quenches terminating at the critical point [Xf = 2, 
Figs.Elc) and Eld)], and except for the largest system sizes, 
display a behavior that is qualitatively similar to the one seen 
in quenches to the extended regime. Namely, they exhibit a 
power-law-like decrease of Am GGE and Arj GGE with increas- 
ing system size. However, a tendency towards saturation can 
also be seen in the differences for the largest system sizes. 
These can be attributed to the failure of the observables to 
relax to stationary values for the times considered here (see 
Figs. [T] and EJ. Hence, as long as relaxation is achieved, the 
GGE provides a good description of observables also at the 
critical point. The GE, on the other hand, fails to describe 
and rj a (as it does in the extended regime). 

The quenches to the localized phase [Xp = 3, Figs.Ele) and 
Elf), and X F = 4, Figs. Etg) and Eth)] exhibit a very differ- 
ent scaling of Am GGE and Atj gge from that observed in those 
to the extended regime and the critical point. One can see 
in the corresponding panels in Fig. [7] that, for Xf = 3 and 
Xf = 4, Af7i GGE and Arj GGE are almost constant with increas- 
ing system size, in the same way (up to an offset) that Am GE 
and At7 ge are. This makes evident that the GGE descrip- 
tion breaks down in the localized phase, in a similar way that 
standard statistical ensembles fail, in general, to describe in- 
tegrable systems after relaxation. We should note, however, 
that the GGE predictions are closer to the long-time averages 
than the ones provided by the GE, as expected given the larger 
number of constraints imposed in the former ensemble. 

We have also studied the scaling of the differences Ati GGE 
and An GE for all parameter regimes depicted in Fig. El We 
find that An GE behaves similarly to Am GE and Arj GE , i.e., 
it saturates to finite values with increasing system size. On 
the contrary, An GGE exhibits a qualitatively different behavior 
from Am GGE and A?7 GGE . Independently of Xf, we find that 
A« GGE is very small and almost size independent. This can 
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FIG. 8: (Color online) Scaling of Arc GGE and An GE with increasing 
system size. Results for An GGE are reported for time averages cal- 
culated using different numbers of time steps. Since for the infinite 
time average An GGE = 0, the number of time steps used in the finite 
average determines the result. The continuous (red) line shows an 
average over 99 steps with T £ [9.9 x 10 5 , 10 6 ], the dashed (green) 
line an average over 990 steps with 16 [9x 10 5 , 10 6 ], and the dot- 
ted (blue) line an average over 9900 steps with T £ [10 4 , 10 6 ]. The 
dash-dotted (red) line shows A;; GE for an average over 9900 steps. 

be understood because is a property that is shared by hard- 
core bosons and noninteracting fermions, and, by construc- 
tion, the infinite time average of one-body fermionic observ- 
ables is given by the GGE. Since for the infinite time average 
An GGE = 0, this quantity is strongly affected by the width of 
the time interval used to calculate the time averages as well 
as by the number of time steps used. Evidence of this de- 
pendence is presented in Fig.|8]for quenches with Xf = 1 and 
Xf = 4 (the results for other values of Xf are qualitatively 
similar). 



One-particle correlations 

The three observables we have studied throughout this work 
provide complementary information about one-particle corre- 
lations. Two of those observables («,■ and mj) are currently ac- 
cessible in ultracold-gas experiments. In order to conclude our 
study, and to make contact with the discussion in Ref. we 
also directly analyze the behavior of one-particle correlations. 
Note that py is a complex Hermitian matrix, and this is why 
rij, m^, and r\ a are all real quantities. 

In Fig. [9] we show how the absolute value of py decays 
when i is fixed to be the central site in the lattice (i = L/2) and 
j moves towards the boundaries. Results are presented for 
two different initial states for quenches towards the extended, 
critical, and localized regimes, for different times (as well as 
for the time average), and within the GGE and the GE. The 
behavior of py in the initial state (in equilibrium py is real) 
reflects the nature of the ground state in the extended and lo- 
calized phases. In the former, one-particle correlations exhibit 
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FIG. 9: (Color online) Decay of the absolute value of py for i = 500 
and j > 500 in a system with L = 1000. The time average was taken 
over 437 steps with % 6 [10 4 ,5.37 x 10 4 ]. The results depicted are the 
absolute values after taking those time averages [py(t) is complex]. 



a power-law decay (fin °= 1/ — no matter the value of 
X, while in the latter they decay exponentially ll52ll . 

The quenches towards the extended phase [Figs. |9{a) and 
[9jb)] exhibit clear similarities no matter the value of Xj. We 
find the following: (i) |py | is very similar, but not the same, for 
T = 100, T = 1000, and the time average, (ii) The time average 
and the GGE results show an excellent agreement with each 
other, (iii) py exhibits a faster, and featureless, exponential 
decay in the GE. This is all consistent with our conclusion 
that the GGE provides an adequate description of one-particle 
observables after relaxation, while the GE fails to do so in this 
regime. 

Figures|9jc) and|9ld) depict results for quenches to the crit- 
ical point. In this case, due to the slow relaxation dynamics 
discussed before, the values of |py| at different times differ 
from each other and from the time average. The time-averaged 
results can be seen to be closest to the GGE prediction and are 
clearly distinct from those in the GE. Calculating the time av- 
erages for later times (not depicted) does improve the agree- 
ment between those averages and the GGE predictions, re- 
vealing a picture similar to the one obtained for quenches to 
the extended phase in Figs.[9|a) and|9jb). 

Results for quenches to the localized phase are presented in 
Figs.[9|e) and|9|f). Once again, |py| at different times differ 
from each other and from the time average. The latter is also 
different (although quite close for the quench A/ = 8 —> Xf = 
4) from the GGE predictions. This is compatible with our pre- 
vious findings that the GGE fails to describe and rj a after 
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relaxation in this regime. Further understanding of the behav- 
ior seen for these quenches can be gained by analyzing the 
case in which Xf °°, so that Hamiltonian (Q~|) can be writ- 
ten as H = £jfij, where Ej is the local chemical potential in 
each site. It then follows that 

Py(T) = (W(T)\b}bj\V(>u)) « Ptj(0yW, (12) 

which means that if one quenches deep inside the localized 
phase, |py(t)| f« Py(0), i.e., correlations present in the initial 
state are preserved, similarly to what we see in Fig.|9je). 

We note that our results in Figs. |9je) and|9jf) are similar to 
the ones reported in Fig. 3 in Ref. 13711 for two-point correla- 
tions of the order parameter. However, the contrasts between 
Figs.|9je) and|9jf) and Figs. [9la) and|9jb) make evident that 
the failure of the GGE in disordered systems is a consequence 
of localization and not of the breaking of translational symme- 
try. Our results also make clear the importance of computing 
time averages, for complex quantities such as py, before com- 
paring with the predictions of the GGE description. 

V. SUMMARY 

In this work, we studied the dynamics and description af- 
ter relaxation of hard-core bosons in one-dimensional lattices 
after a sudden change of the strength of an additional quasi- 
periodic potential. This model features two distinct regimes, 
an extended regime for weak quasi-periodic potentials and 
a localized regime for strong quasi-periodic potentials. Our 
analysis has shown that the approach of observables towards 
their time-independent values after relaxation comes close to 
following a power law. For the finite system sizes studied, 
all observables reach their time-independent values within 
the considered time scales. The sole exceptions were the 
quenches towards the critical point, where the dynamics was 
found to be slower and time-independent values of the ob- 
servables were not reached for the largest lattices. We have 
argued that, in most of the cases analyzed, the times required 
for the observables to reach their stationary values increase as 
a power law with the system size. 



We further compared the long-time average of observables 
with statistical descriptions provided by the GE and the GGE. 
The GE failed to describe all observables after relaxation in 
the quenches considered, as expected since these systems are 
integrable. The GGE, on the other hand, was found to pro- 
vide a good description of observables after relaxation in the 
extended phase, and at the critical point, whenever observ- 
ables became time independent (up to vanishingly small fluc- 
tuations). The scaling behavior in these two cases suggests 
that, in the thermodynamic limit, the GGE results are identi- 
cal to those after relaxation. On the contrary, in the localized 
regime, we have found that the GGE fails to describe observ- 
ables that depend on nonlocal correlations (such as and 
T]a) after relaxation, and that this picture does not change with 
changing system size. The time average of the density, on the 
other hand, was shown to be well described by the GGE in all 
regimes. 

From the outcome of this study, as well as from the results 
in Refs. |[H-i3l [l7ll . we conclude that localization, and not the 
breaking of translational symmetry as proposed in Ref. ll3~7ll , 
can lead to the breakdown of the GGE description. Our work 
also poses the question of whether modifying the GGE by us- 
ing a different set of conserved quantities (here we used the 
occupation of the single particle eigenstates of the noninter- 
acting fermionic system to which hard-core bosons can be 
mapped), or adding further conserved quantities, would allow 
one to describe time averages of observables in the localized 
regime. Recent work on finding optimal sets of conserved 
quantities may shed light on these questions ll58ll . 



Acknowledgments 

This work was supported by NSF under Grant No. OCI- 
0904597 and by the Office of Naval Research. We thank E. 
Khatami, K. He and T. Wright for helpful discussions and 
comments on the manuscript. 



[1] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. 

Lett. 98, 050405 (2007). 
[2] M. Rigol, A. Muramatsu, and M. Olshanii, Phys. Rev. A 74, 

053616 (2006). 

[3] A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. Lett. 106, 

140405 (2011). 
[4] M. A. Cazalilla, Phys. Rev. Lett. 97, 156403 (2006). 
[5] P. Calabrese and J. Cardy, J. Stat. Mech. p. P06008 (2007). 
[6] M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. 

Rev. Lett. 100, 030602 (2008). 
[7] T. Barthel and U. Schollwock, Phys. Rev. Lett. 100, 100601 

(2008). 

[8] M. Eckstein and M. Kollar, Phys. Rev. Lett. 100, 120404 
(2008). 

[9] M. Kollar and M. Eckstein, Phys. Rev. A 78, 013626 (2008). 



[10] A. Iucci and M. A. Cazalilla, Phys. Rev. A 80, 063619 (2009). 
[11] D. Rossini, A. Silva, G. Mussardo, and G. E. Santoro, Phys. 

Rev. Lett. 102, 127204 (2009). 
[12] D. Rossini, S. Suzuki, G. Mussardo, G. E. Santoro, and 

A. Silva, Phys. Rev. B 82, 144302 (2010). 
[13] D. Fioretto and G. Mussardo, New J. Phys. 12, 055015 (2010). 
[14] A. Iucci and M. A. Cazalilla, New J. Phys. 12, 055019 (2010). 
[15] J. Mossel and J.-S. Caux, New J. Phys. 12, 055028 (2010). 
[16] P. Calabrese, F. H. L. Essler, and M. Fagotti, Phys. Rev. Lett. 

106, 227203 (2011). 
[17] M. A. Cazalilla, A. Iucci, and M.-C. Chung, Phys. Rev. E 85, 

011133 (2012). 

[18] M. Rigol and M. Fitzpatrick, Phys. Rev. A 84, 033640 (201 1). 
[19] M.-C. Chung, A. Iucci, and M. A. Cazalilla, New J. Phys. 14, 
075013 (2012). 



12 



[20] M. Kormos, A. Shashi, Y.-Z. Chou, and A. Imambekov, 

arXiv: 1204.3889 (2012). 
[21] K. He and M. Rigol, Phys. Rev. A 85, 063609 (2012). 
[22] R Calabrese, R H. L. Essler, and M. Fagotti, J. Stat. Mech. p. 

P07016 (2012). 

[23] P. Calabrese, R H. L. Essler, and M. Fagotti, J. Stat. Mech. p. 
P07022 (2012). 

[24] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and 

M. Rigol, Rev. Mod. Phys. 83, 1405 (2011). 
[25] T. Kinoshita, T. Wenger, and D. Weiss, Nature (London) 440, 

900 (2006). 

[26] S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, 
U. Schollwock, J. Eisert, and I. Bloch, Nature Phys. 8, 325 
(2012). 

[27] E. T. Jaynes, Phys. Rev. 106, 620 (1957). 
[28] E. T. Jaynes, Phys. Rev. 108, 171 (1957). 
[29] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991). 
[30] M. Srednicki, Phys. Rev. E 50, 888 (1994). 
[31] M. Rigol, V. Dunjko, and M. Olshanii, Nature (London) 452, 
854 (2008). 

[32] M. Rigol and M. Srednicki, Phys. Rev. Lett. 108, 110601 
(2012). 

[33] J. Berges, S. Borsanyi, and C. Wetterich, Phys. Rev. Lett. 93, 
142002 (2004). 

[34] M. Moeckel and S. Kehrein, Phys. Rev. Lett. 100, 175702 

(2008) . 

[35] M. Moeckel and S. Kehrein, Ann. Phys. (N.Y.) 324, 2146 

(2009) . 

[36] M. Kollar, F. Wolf, and M. Eckstein, Phys. Rev. B 84, 054304 
(2011). 

[37] T. Caneva, E. Canovi, D. Rossini, G. E. Santoro, and A. Silva, 



J. Stat. Mech. p. P07015 (2011). 
[38] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Ann. Phys. 321, 
1126 (2006). 

[39] V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155 1 1 1 (2007). 

[40] A. Pal and D. A. Huse, Phys. Rev. B 82, 17441 1 (2010). 

[41] E. Khatami, M. Rigol, A. Relano, and A. M. Garcfa-Garcia, 

Phys. Rev. E 85, 050102(R) (2012). 
[42] C. Gogolin, M. P. Miiller, and J. Eisert, Phys. Rev. Lett. 106, 

040401 (2011). 

[43] S. Aubry and G. Andre, Ann. Israel Phys. Soc. 3, 133 (1980). 
[44] J. B. Sokoloff, Phys. Rep. 126, 189 (1985). 
[45] T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940). 
[46] E. Lieb, T. Shultz, and D. Mattis, Ann. Phys. (NY) 16, 406 
(1961). 

[47] P. Jordan and E. Wigner, Z. Phys. 47, 631 (1928). 

[48] M. Rigol and A. Muramatsu, Phys. Rev. Lett. 94, 240403 

(2005) . 

[49] O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956). 

[50] A. J. Leggett, Rev. Mod. Phys. 73, 307 (2001). 

[51] A. M. Rey, 1. 1. Satija, and C. W. Clark, Phys. Rev. A 73, 063610 

(2006) . 

[52] K. He, 1. 1. Satija, C. W. Clark, A. M. Rey, and M. Rigol, Phys. 

Rev. A 85, 013617 (2012). 
[53] N. Nessi and A. Iucci, Phys. Rev. A 84, 063614 (201 1). 
[54] M. Rigol and A. Muramatsu, Mod. Phys. Lett. 19, 861 (2005). 
[55] M. Rigol, Phys. Rev. A 72, 063607 (2005). 
[56] P. Harper, Proc. Phys. Soc. London A 68, 874 (1955). 
[57] L. C. Venuti and P. Zanardi, arXiv: 1208. 1121 (2012). 
[58] M. Olshanii, arXiv: 1208.0582 (2012). 



